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A METHOD FOR COMPUTING CHEMICAL -EQUILIBRIUM COMPOSITIONS 
OF REACTING-GAS MIXTURES BY REDUCTION 
TO A SINGLE ITERATION EQUATION 

By Wayne D. Erickson, Jane T. Kemper, 
and Dennis O. Allison 
Langley Research Center 

SUMMARY 

A method for computing the equilibrium chemical composition of a reacting gas 
mixture is presented wherein all the equilibrium and mass-balance equations are com- 
bined algebraically into one equation involving a single imknown. A comparison between 
this method and a widely used more general method shows that the unit computational 
speed of the present method can be as much as two orders of magnitude faster. This 
great advantage in computational speed is achieved at the expense of generality. The 
present method for computing equilibrium chemical compositions should be useful 
when applied to systems of intermediate complexity and where a large number of com- 
putations are required. 


INTRODUCTION 

A number of methods have been developed for computing the equilibrium chemical 
composition and thermodynamic properties of reacting gas mixtures. These methods 
employ various iteration schemes. Although the machine computation time required to 
obtain a single solution may be small, the machine computation time in flow -field prob- 
lems that require a very large number of equilibrium composition calculations may be 
excessive. In that case, a series of empirical curve fits of the thermodynamic prop- 
erties is sometimes used. 

A study of three of the most widely used general methods for computing chemical- 
equilibrium compositions was carried out by Zeleznik and Gordon (ref. 1.) The study 
considered the methods of Brinkley (ref. 2), Huff, Gordon, and Morrell (ref. 3), and White, 
Johnson, and Dantzig (ref. 4) with regard to any computational advantage of one method 
over the other methods. The comparison showed that none of the methods indicated any 
significant computational advantage over the other two. 



The present work outlines a method for computing the equilibrium chemical com- 
position wherein the complete set of equilibrium and mass-balance equations are reduced 
to one equation involving a single unknown. A comparison between the present method 
and the more general method of White will show that a great advantage in computational 
speed is realized in the present method. As will be shown, this advantage is achieved at 
the expense of generality. 


SYMBOLS 


a^ y. mass of kth element per mass of ith species 

aj^ number of atoms of kth element per molecule of ith species 


h. 

RT 

\,k+l 


c^(n) 




RT 



nondimensional Helmholtz free energy of ith species 

. .Xj^ nondimensional Helmholtz free energy of reacting gas mixture 

mass ratio of kth element to the (k + l)th element in the mixture, regard- 
less of chemical species in which they appear 

total number of atomic weight units of element k per mass of mixture 

mass fraction of species i in mixture 

nth approximation of Cj^ 

indexed coefficients in equation (6) 

function defined by equation (6) 

derivative of (see eq. (8)) 

nondimensional Gibbs free energy of ith species 

. .XjN nondimensional Gibbs free energy of reacting gas mixture 
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lA 


RT 


■S=,i 
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N 

P 

R 

T 

Wi 

W 


X 


X. 


|8 


^out 

*T 


nondimensional standard -state free energy of formation of species i 
at temperature T 

integers 

equilibrium constant in terms of partial pressures for reaction j 
unit computational time, time/case 
number of cases computed per run 
pressure 

universal gas constant 
temperature 

molecular weight of ith species 

molecular weight of mixture 

mole number of species i, Moles of species i 

Mass of mixture 

mole number of mixture, Total moles of mixture 

Mass of mixture 

chemical symbol of species i where, for example, in equation (4) 

Xj = H, Xg = O, Xg = Hg, etc. 

mass-balance quantity defined by equation (22) 

actual computational time (this time is exclusive of input and output times) 

input time 

output time 

total machine time 
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V- ■ stoichiometric coefficient of reactant species i in ith reaction 

v\ . stoichiometric coefficient of product species i in jth reaction 

p density of mixture 

Tj modified equilibrium constant defined by equation (5) 

Subscripts: 

i ith species 

I total number of species in mixture 

j jth reaction 

J total number of independent reactions that can be written for mixture 

k kth element 

K total number of elements in mixture 

I power series index in equation (6) 

L highest order power in equation (6) 

EQUILIBRIUM COMPOSITION EQUATIONS 

This section is concerned with the development of a technique for computing the 
equilibrium chemical composition of a reacting mixture at a given temperature, density, 
and set of mass restraints, wherein the solution for a single equation involving only one 
composition term is obtained. The solution for all the remaining species then are found 
by direct substitution. The general equations for an equilibrium reacting gas system 
involving an arbitrary number of species and chemical elements are presented, after 
which a specific set of equations are written for a reacting mixture of hydrogen and air. 

General Equations 

In an equilibrium reacting gas mixture containing I chemical species, which in 
turn are composed of K chemical elements, there are J independent and simultaneous 
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chemical reactions that can be written, such that I = J + K. Application of the conser- 
vation of mass for each of the kth chemical elements leads to the following K equations: 


and 


I 

1 

Y =^k,k+l (k=l,2,. . .K-1) (1) 

Z^k+l'^i 

i=l 


I 



i=l 


( 2 ) 


The quantity a^^ represents the total mass of atoms of the kth element per mass of a 
molecule of the ith chemical species, a^ is the corresponding quantity for the 
(k+l)th element, and is defined as the mass ratio of the kth element to the 

(k+l)th element in the mixture, regardless of the chemical species in which the element 
appears. 


There are also J equilibrium expressions which can be written as 


I 


¥ 


T 




(j — 1,2,. . .J) 


(3) 


which correspond to the J independent chemical reactions containing species 


I 

1=1 







a . 1,2,. . .J) (4) 


The quantity Tj in equation (3) is a modified equilibrium constant which is a function 
of T and p. 


"j"' = i ■'i.r "1,J IT wr*’ 

i=l 


.-V. . 

] 1,] 


(j = 1,2,. . .J) 


(5) 


where K„ . 

P,3 


is the equilibrium constant in terms of partial pressures for the jth reaction. 
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Equations (1), (2), and (3) constitute I independent equations involving as many 
vuiknown Cj^. The present technique for computing the chemical composition at a speci- 
fied temperature, density, and mass restraints (specific mass restraints are reflected 
through j^, a. and bj^ of equation (1)) first entails writing down the 
K equations suggested by equations (1) and (2). These K equations contain all 
I unknowns, namely, the mass fractions Cj. The J equilibrium equations given by 
equation (3) are then substituted into these K equations for the mass constraints in 
order to reduce the I unknowns by J so that a modified set of K equations is 
obtained containing only K unknowns, for example, CpC 2 ,. . These resultir^ 

K equations are then combined algebraically in such a way as to obtain a single equation 
containing only one unknown, for example, Cj^, where Cj^ is the mass fraction of an 
arbitrarily chosen species 1. (See appendix for the details of the algebraic reduction.) 
This single equation can be expressed as 

L 

*'(<=!) = ^ Oj'!* = « <6) 

1=0 

where the coefficients d^ are fimctions of a^ j^, a^ bj^ y i/! j, and Tj 

but are independent of the mass fractions c.. It follows that the coefficients d^ can be 
computed once the temperatvire, density, and set of mass constraints are imposed. The 
calculated values of d^ can then be used in equation (6) along with the application of 
Newton's iteration scheme to find Cj^. 

Newton's iteration scheme for determining Cj^ can be written as 

F(cj(n)) 

c^(n + 1) = Cj(n) - ^ (7) 

^ 1 F’(c^(n)) 

where Cj^(n) is the nth approximation of Cj^ and the derivative of F^Cj^J is 

L 

F'(Ci) = 2 idjC,'-' (8) 

1=1 


The iteration scheme indicated by equation (7) is repeated until the required accu- 
racy for Cj is achieved. Equation (6) indicates that a number of roots are possible 
for c J. However, the correct value of Cj^ must be within the rai^e 0 s ^ 1. If a 
root outside this range is obtained through the iteration process, it is eliminated and a 
different initial choice of Cj^ is made to restart the iteration process. After a solution 
for Cj^ is found within the range 0 i c^ = 1 through equation (7), the mass fractions 
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for all the other species, found by direct substitution into the inter- 

mediate equations that lead to equation (6). (These intermediate equations are not shown 
for the general case but will be indicated for the specific example for a hydrogen-air 
model which follows.) Each value of the mass fractions, 02 , 0 ^,. . .Cj, is examined to 
make certain it falls in the range from 0 to 1. The sum of all the mass fractions must 
also be unity. There is only one set of values for c^ which satisfies these conditions. 


Specific Equations for a Hydrogen-Air Gas Model 

The application of the present technique for computing the chemical composition is 
best illustrated by a specific example of intermediate complexity. Consider an equilib- 
rium system resulting from the combustion of pure hydrogen with air. Assume that air 
is a 4 to 1 molar mixture of N2 and O2 wherein the N2 is taken as inert. The following 
seven species (I = 7), H, O, H2, O2, OH, H2O, and N2 with species index numbers 
i = 1,2,. . .7, respectively, are assumed to appear in the system; the molecular weight of 
each of these species is taken to be 1, 16, 2, 32, 17, 18, and 28, respectively. These 
species are composed of three elements (K = 3), H, O, and N with element index numbers 
k = 1, 2, and 3, respectively. The two specific equations that result from equation (1) 
with k = 1 and k = 2 are, therefore, 


Cl + Co + — Cc + - c^ 

1 3 17 5 9 6 

bi 2 

16 8 
Co + c. + — Cc + - 

2 4 17 5 9 6 


and 


16 8 

Co + c - + --- c,- + - c^ 

2 4 17 5 9 6 


= b 


2,3 


2 

7 


( 9 ) 


( 10 ) 


where b^^ 2 the arbitrary mass ratio of hydrogen to oxygen in the system and 
2 

^2 3 ” 7 from the assumption that air is a 4 to 1 molar ratio of N2 to O2. In 

addition to these two equations, equation (2) applied to the present example leads to 


Cf + C2 + C3 + C4 + Cg + C6 + C7 = 1 (11) 

In addition to these three equations, there are J = I - K = 4 independent equilib- 
rium expressions which result from the following independent chemical reactions. 


7 



(3 = 1) 


Hg = 2H 

Og =20 (j = 2) 


OH = O + H (j = 3) 

HgO = 2H + O (j = 4) 

Equation (3) then leads to the following equations where j has been set equal to 1, 2, 

3, and 4, 


3 " "^l^^l^ 

(12) 

4 '^2*^2^ 

(13) 

= nxgCjCg 

(14) 

= 9t^Cj^C2 

(15) 


(The arbitrary inclusion of the factors 17 and 9 in equations (14) and (15) makes for a 

somewhat simpler set of values for d^, as will be shown. The only effect is to modify 

the relation between t. and K . shown by equation (5).) The quantities t. in these 

J P> J 3 

equations are related to the equilibrium constants in terms of partial pressures Kp j 

for each reaction through equation (5) and can be expressed explicitly as * 


2pRT 

■■S,! 

(16) 

pRT 

(17) 

pRT 

(18) 

W 

00 

INO 

(19) 


where Kp j are functions of temperature alone. (The factors 17 and 9 in equations (14) 
and (15) have been accoimted for in writing down equations (18) and (19). If these factors 
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were not used, then in what follows, it would be necessary to replace Xg by and 

T4 by 9x4.) 

It follows then that for a given temperature, density, and mass ratio of hydrogen to 
oxygen b^ g? equations (9) to (15) constitute seven independent equations involving seven 
unknown quantities, 04,02,. . .Crj. Equations (10) and (11) are combined to eliminate Crj 
and equations (12) to (15) are then employed to eliminate Cg, C4, Cg, and Cg to give 

9x 2C2^ + (74x404^ + 146xgC4 + 9^C2 + 2x404^ + 2C4 - 2 = 0 (20) 

Equation (9) can also be reduced to contain only C4 and C2 by substitution of equa- 
tions (12) to (15) so that 

^X4C4 + ^1^1 ^1 ~ ^ ^ (21) 


where 




2b 


lii. 


2b 1,2 + 9 


(22) 


Equation (21) rearranges to give an explicit expression for C2 in terms of C4 only. 


^ - "^^^4 - C4 

X4C12 + XgC4 


(23) 


This explicit expression for C2 in terms of C4 is then substituted into equation (20) 
to obtain a single equation that contains C4 only. With some rearrangement, this 


expression containing C4 only can be written as 

6 

>'{‘=1) = I Vl' = » <2“) 

z=o 

where the coefficients d^ are 

dg = (24a) 

di = -( 2 x 2 - Xg^)3 (24b) 

^2 = (^2 - - (2'"i^ 2 - '^4)^ + - 0^3^ (24c) 
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(24d) 


^3 = (2^1 ^2 - ^rs - - ^4) + f - 1)^3 ^4 

^4 = (^1^^2 - ^1^4 - 16^1^3^ - 24r3T4) + |(37^ - l)r^^ 
dg = -8{3t^t^t^ + 7 ^ 2 ) 


(24e) 

(24f) 

(24g) 


A solution for Cj can be obtained by application of equation (7) (Newton's itera- 
tion scheme) which yields 


6 

X d^c^l(n) 

Cj(n+1) = c, (n) (25) 

W^c^l"^(n) 



wherein Cj^(n) is the nth approximation for Cj and Cj(n+1) is the next better approx- 
imation in this iteration scheme. The number of iterations required is, of course, 
dependent upon the accuracy desired. 


For this specific example in which it has been assumed that air is a 4 to 1 molar 
mixture of N 2 and O 2 , and the mass ratio of hydrogen to oxygen is denoted by bj^ 2 > the 


possible range of c^ is O^Cj^g 


1.2 


H ^9 

Di rt + — ' 

1.2 2 


After a value of Cj^ in this range has been 


obtained through equation (25) with the desired accuracy, the value of C 2 can be com- 
puted directly from equation (23). Equations (12) to (15) are then employed to compute 
Cg to Cg, respectively. The value of Crj is then computed by difference from equa- 
tion (11). It should be noted again that all the values of c^^ must each be within the 
range 0 = c^ = 1. 

The procedure for applying this technique to other systems would be the same as 
indicated in this specific case. As the system becomes more involved with respect to the 
number of species and reactions, the order of the series equivalent to equation (24) 
increases and it becomes increasingly more tedious to determine the coefficients d^. 
This present technique is probably best suited to systems of intermediate complexity. 
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COMPARISON WITH OTHER METHODS 


The purpose of this section is to compare the present method with other existing 
methods for computing the equilibrium chemical composition of reacting gas mixtures 
with particular regard for the speed of computation. As already mentioned, a study by 
Zeleznik and Gordon (ref. 1) has shown that none of the three most widely used general 
methods offer any significant computational advantages over the other. A comparison 
between only one of these methods and the present method is therefore all that is needed. 
The White method (ref. 4) which employs a free -energy minimization technique was 
selected for the purpose of making this comparison. The comparison to be made here is 
limited to determining the chemical composition at a given temperature and density for 
the two methods. It is recognized that solutions to problems are often required when two 
other thermodynamic quantities are specified. Where the latter is desired, the present 
method could become a part of an iteration scheme based on two thermodynamic quanti- 
ties. It is believed that the most time-consuming step in all such schemes is the deter- 
mination of the equilibrium composition; therefore, only the timing of this process will be 
considered. 


Free -Energy Minimization Technique 

The determination of the equilibrium chemical composition through the free -energy 

minimization technique requires the standard free energy of formation AF?. of each 

1,1 

species in the mixture at the temperature of the mixture. The free energy of the mix- 
ture, which is a minimum at the equilibrium chemical composition, can be written as 




fAF 






RT 


+ logg p + logg 



(26) 


where . .Xj represent a set of I mole numbers, Xj^ is defined as the moles of 

species i per gram of mixture, p is the total pressure expressed in atmospheres, and 

I 

X = ^ Xj (27) 

i=l 


The equilibrium chemical composition is determined when a set of mole numbers 
XjjXg,. . .Xj produces a minimum value for the free energy of the mixture 

F 

^;j^Xj,X 2 ,. . .Xj^ given by equation (26) and at the same time satisfies the following mass 
balance expressions: 
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(k=l,2,...K) (28) 


I 

i=l 

where represents the number of atoms of element k contained in a molecule of 

species i, and bj^ is a constant equal to the total number of atomic weight units of 
element k per gram of mixture. Equation (26) when minimized within the limitations 
of equation (28) yields the chemical composition at a fixed T and p. Since the present 
method was developed for a fixed T and p, and a comparison between the two methods 
is desired, equation (26) is rewritten in terms of T and p. The equation of state of the 
mixture can be expressed as 


P = P W T = pxRT 

Substitution of equation (29) for the pressure into equation (26) yields 

I 

i(x,,X2,. . .Xi) = 1 Xj 
i=l 


AFfOi 

+ logg p + logg(x. ) + logg RT 


(29) 


(30) 


wherein the units of R, where it is not a part of a nondimensional quantity, must be such 
that pressure takes on units of atmospheres. However, the appropriate quantity to min- 
imize for a fixed T and p is no longer ~^Xj,X 2 ,. . .Xj\, the Gibbs free energy, but 

rather ’ *^l)’ Helmholtz free energy. The relationship between the 

^i 

Gibbs molar free energy of species i, — and the Helmholtz molar free energy of 

. , A. , 

species 1 , — IS 



h. 

RT 


•+ 1 


(31) 


Therefore, the expression to be minimized to find the equilibrium composition for a 
fixed T and p becomes 


A 

RT 



AF 


^i 


f.i 


+ logg p + logg(Xj^ + logg RT 


1 


(32) 


Since a machine program based on the minimization of equation (26) at constant p 
and T was available from a previous study, a program was not developed to use equa- 
tion (32) at constant p and T directly. Instead, equation (26) was replaced by 
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equation (30) in the existing program and the input included fixed values of p and T. 
Since the existing program was set up for constant p, it follows that each iteration step 
is carried out for a constant pressure equal to that given by equation (29), but this pres- 
sure is modified from one iteration step to the next in direct proportion to x. As the 
result approaches the correct solution, the pressure stabilizes. The computation time 
required for this solution was found to be essentially the same as that required when a 
given temperature and pressure were used directly in the program that utilized equa- 
tion (26). The numerical results were also the same. 


Input Quantities for Comparison 

The input quantities for both the present method and the free -energy minimization 
technique include T, p, and a set of values for the molar standard -state free energy of 


formation of each species 

AF,"' 


af: 


at T. The free -energy method makes use of the quan- 
RT 

AF° 


f i f i 

tity directly in equation (30) whereas the present method utilizes 2 - in com- 

RT RT 

puting equilibrium constants Kp j which are required in equation (5) for computing a set 


of values of Ty 


AF 


The relation between K . and 

P,1 


lA 


RT 


IS 


Kp,j = exp 


4 I ‘>V RT 

i=l 


(j = 1,2,. . .J) (33) 


In addition to these common quantities, the direct method requires an initial approxima- 
tion for the mass fraction of H, that is, Cj^; whereas the free -energy method requires 
an initial approximation for a set of mole numbers of all species. This set of mole num- 
bers must, of course, be chosen to be consistent with the mass balance expressions 
given by equation (28). 


The temperature and density chosen for comparison are T = 4000° K and 

-5 n i 

p = 3.0 X 10 g/cc. At 4000° K the free -energy quantities from reference 5 are 

RT 


Species 


RT 

Species 

RT 

0 

H 

-0.46052 

0 

OH 

-.52269 

0 

-.39144 

H 2 O 

-.54802 
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A stoichiometric mixture of H 2 and air, wherein air is assumed to be a 4 to 1 molar 
mixture of Ng to O 2 , is used to define the mass balance. This definition requires that 
the quantity b^^ 2 equation (22) be set equal to 1/8 for the present method, and for the 
free-energy method, a suitable initial set of mole numbers is xS- = 0.0270270, 

INfl 


Xq = 0.0067568, x^ = 0.0135135, and Xq = x^ = Xqjj = q = 10“^. This set corre- 

spends to a 1;2:4 molar mixture of 02.*H2:N2 in a reference state. Because the program 
for the free-energy minimization technique requires initial nonzero values for all species 
considered in the mixture, the reference values for the mole numbers of O, H, OH, and 
H 2 O were chosen to be outside the eight significant number range but nonzero. The 
present method requires an initial estimate of the mass fraction of only one species and 
a value of zero is chosen for the mass fraction of atomic hydrogen, Cj^. The results of 
the calculations by the two methods are discussed in the next section. 


RESULTS AND DISCUSSION 


The present method and the free-energy minimization technique for computing the 
chemical composition of the aforementioned hydrogen-air mixture were programed in 
FORTRAN IV. A number of runs were made using these two programs on an IBM 7094 
computer at the Langley Research Center. These computations were made at a temper- 
ature of 4000° K, total density of 3 x 10"^ g/cc, and a stoichiometric mixture of hydrogen 
in air ^4 to 1 molar ratio of N 2 to O 2 ). The accuracy criterion for each method was 
adjusted so that the hydrogen atom mass fraction would be calctilated to five -place accu- 
racy for both methods. The calculated results for these two methods, expressed in terms 
of mass fraction of each species, are listed in table I. The variations in the accuracy 
between the two methods are indicated by underscoring. 


TABLE I.- EQUILIBRIUM MASS FRACTIONS FOR SPECIES IN STOICHIOMETRIC. 
HYDROGEN -AIR MIXTURE FOR TWO COMPUTATIONAL METHODS 
|t = 4000° K and p = 3 x g/c^ 


Species 

■■ 

Present method 

Free-energy minimization technique 

No 

0.75675675 

0.75675599 

2 

.16356338 x 10~^ 

.16356545 X 10"^ 


.35859237 X 10~^ 

.35859171 X 10"^ 

0 

.17050637 

.17050748 

H 

.213864;^ X 10"^ 

.21386453 x 10"^ 

OH 

.27447568 X 10"^ 

.27447716 X 10"^ 

HgO 

.39606028 x 10"^ 

.39606200 X 10"^ 
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After it had been established that the two methods give the same numerical result 
within the same order of accuracy, a number of computational runs were made to com- 
pare the speed of calculation for each method. The timing technique available for this 
study permitted only a determination of the total running time for an arbitrary number of 
cases. This total time is composed of three parts: the input time the actual com- 
putational time t^, and the output time The computational time for N identical 

cases can be expressed as the product of N and the computation speed m 

tg = Nm (34) 


and the total machine time is 


^T = *in + 


Nm + t 


out 


(35) 


The quantity of interest in this study is m, the actual computational speed for each 
method; however, only trp is measured directly while N is an input value. One way 
of determining m is to obtain values of t,p over a range of N for conditions such 
that Nm » t^^^ and Nm » In this study the program was allowed to run for a 
given N without an output so that t^^^ = 0. This procedure was followed after the two 
methods were shown to yield the same solution to the same order of accuracy as indicated 
in table I. The second condition, Nm » was achieved by measuring trj, for a very 
large number of cases, N. 

The results for a series of timed machine calculations for both the present method 
and the free -energy minimization technique are presented in figure 1. The total machine 
running time t,p is plotted as a function of the number of cases run N for the free- 
energy minimization technique and for the present method. The previously mentioned 
initial set of mole numbers were used for the free -energy minimization technique which, 
in turn, required 12 iterations to find the solution. These results are represented by the 
uppermost curve where N ranges from 500 to 8000. It is seen that the slope of this 
curve has reached unity for N = 6000; thus, in this region, the unit computational time is 

m = ^ (Nm » = 0) (36) 

The unit computational time m is, of course, independent of both t,p and N so that 
the value of m determined in this manner is the time per unit case for any number of 
cases run. The value of m for the free -energy minimization technique requiring 
12 iterations is seen from figure 1 to be 252 msec/case. The deviation from unit slope 
for the lower values of N on this curve is due to the increasing importance of the input 
time t^^ relative to Nm as N decreases. 
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Comparison of running times for free-energy minimization technique and present method. 



The free -energy minimization technique was also timed for two additional runs for 
which the initial set of mole numbers were chosen to be much closer to the ultimate 
values. The number of cases run N for each of these two runs were 2000. One of these 
runs resulted in 5 iterations, while a run for a still more accurate initial set of mole 
numbers resulted in 4 iterations. These two runs are also shown. All the effect of t.^^ 
has not been eliminated in these two runs so that the value of m for each was estimated 
by drawing in a line of unit slope displaced slightly below each data point to approximate 
the influence of t^^^. The run requiring 5 iterations resulted in a unit computational 
time of 156 msec/case and the run requiring 4 iterations gave m = 108. It should be 
mentioned that the runs that required only 4 or 5 iterations resulted from more accurate 
initial sets of mole numbers. The series of runs in which 12 iterations were required 
is probably a more practical basis for comparison. 

The present method of calculation was also timed for a series of runs for the same 
computational conditions, that is, no output time = 0) a^nd over a range of N with 

N reaching large values. A series of runs were made for each of two different iteration 
schemes. The results for a halving-mode* iteration scheme is shown for the range 
500 i N ^ 8000. The scattering of the data points is due to uncontrollable but slight vari- 
ations in tjj^ from run to run. The total running time even at N = 8000 is still 
strongly influenced by the input time; therefore an accurate evaluation of m is not 
attempted. The value of m would, however, be less than 12 msec/case. An accurate 
value of m could be obtained by making computations for still larger values of N. 

This was not done because the halving-mode scheme is not one of the most timewise 
efficient iteration schemes. Instead, the Newton iteration scheme indicated by equa- 
tion (7) was used to complete the comparison. 

The lowest curve in figure 1 shows the results for the present method using the 
Newton iteration scheme. A series of runs were carried out for an initial value of Cj^ 
equal to 0.001 and another series for an initial value of zero. Within the scatter of the 
data, which again is due to slight variations in tj^^ from rvm to run, the two initial values 


*The halving -mode scheme makes use of equation (24) directly. Initial values of 
Cj equal to 0, 0.5, and 1.0 are used to evaluate The interval in which a sign 

change in noted is then assmned to contain the desired root. The midpoint of 

this new and more limited interval is then computed and used to seek a more narrow 
interval by observing a sign change in F^Cj^J. This procedure is repeated imtil the 
desired accuracy is achieved. If more than one root exists in the region 0 § ^ 1.0, 

the iteration scheme must accoimt for this. There is, of course, only one root c^ that 
satisfies all the other equations whereas all values of Cj are real and positive. 



of Cj produced the same curve. In these series of runs it was necessary to make N 
as large as 10^ in order to satisfy the condition that Nm » tj^^. A slope of unity is 
essentially achieved when N = 10^ so that m = 2.52 msec/case. This unit computation 
time m for the present method using the Newton iteration scheme is two orders of mag- 
nitude less than that for the free -energy minimization technique with 12 iterations. This 
comparison is for the same order of accuracy in the composition for each method, as 
well as the same input quantities. 

The free -energy minimization technique and the other general methods for com- 
puting the equilibrium chemical composition of reacting mixtures can of course be used 
for a large range of reacting gas mixtures with a very large number of species and ele- 
ments. This generality offers a great advantage whereas the present method is limited 
to the species and elements (not their relative amovints) that have been chosen to develop 
the iteration equations. The advantage of the present method over the more general 
methods is realized when machine computation time is important, for example, when a 
finite -difference calculation network requires a very large number of equilibrium com- 
position calculations. The computations necessary to determine the equilibrium chemical 
composition from point to point in certain flow field calculations often require excessive 
machine time so that a series of empirical curve fits of the equilibrium thermodynamic 
quantities are used to achieve a more reasonable machine time. This representation of 
the thermodynamic quantities is limited to the specific gas model for which it was devel- 
oped. If the present method is used, account could be made for the point to point change 
in the elemental distribution due to diffusion by allowing the quantities bj^ 2>^2 3’" ' ‘ 
in equation (1) to be variables determined by the appropriate diffusion equations. 

It is noted that the present method offers considerable advantage over the more 
general method with regard to machine computation time. This advantage has been 
achieved at the expense of loss of generality. The comparison herein was based on a 
hydrogen-air gas model which represents a system of intermediate complexity. Other 
systems of similar complexity could likewise be formulated and used in conjimction with 
the present method. It is possible to apply this technique to more complex systems. It 
does not seem practical, however, to extend this present method to very complex systems 
because the algebraic manipulation becomes overwhelmingly tedious as the number of 
species and/or elements in the system become very large. For this reason, the practical 
application of the present method is probably limited to systems of intermediate com- 
plexity like the system studied or perhaps only somewhat more complex. 
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CONCLUDING REMARKS 


A method for computing the chemical composition of reacting gas mixtures has 
been presented that shows considerable advantage with respect to machine computation 
time. For the particular comparison considered herein, the unit computational speed for 
this method is two orders of magnitude faster than the more general and most commonly 
used computation schemes. This advantage is achieved by reducing the system of equa- 
tions down to one equation containing a single unknown. This reduction to one equation 
involves much tedious algebra, but the procedure is believed to be worthwhile when com- 
putational time must be kept short and when a large number of computations are required 
for the same chemical system. The application of the present method to machine com- 
putations requires that a new set of equations be developed for each chemical system that 
contains a different set of chemical species. The relative total mass of each element in 
the system can be treated as variable input. 

The great advantage in computational speed noted was based on a comparison 
wherein the chemical system was of intermediate complexity. It is believed that chem- 
ical systems of greater complexity would also show this advantage in computational 
speed, but the effort required to reduce the system of equations to the necessary single 
expression may be the limiting factor for very complex systems. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., March 16, 1966. 
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APPENDIX 


REDUCTION OF WORKING EQUATIONS TO A SINGLE EQUATION 
CONTAINING A SINGLE UNKNOWN 


The substitution of the J equilibrium expressions indicated by equation (3) into 
the K conservation of mass expressions indicated by equations (1) and (2) yields 
K equations which contain only K unknown mass -fraction variables. These K equa- 
tions can be written in general form as 

I (n=l,2„ ..K) (Al) 

m=0 


where denotes the mass fraction of a particular species which is one of the K 
mass-fraction variables in these K equations, is the highest order of c^ in the 

nth equation, and A ^ is a polynomial function of the remaining K - 1 mass- 
fraction variables from which c^ is excluded, as well as a function of the known quan- 
titles T., a. j^, I.: ., and 


Equation (Al) can be divided through by either A^ or 

X ^ to give either 

A 

(n=l,2,. ..K) (A2) 


1 + y 

nfel ^0,n 




or 


M -1 

M_ A A 


""n V 

" Z 


m=0 ^n>*^ 


M-c 0 
a 


(n=l,2,. ..K) (A3) 


Consider equation (A2) first with n = 1 and then with n = 2 so that 


Ml 

E Aoi 

m=l 


m J, c = 0 

a. 


(A4) 


and 


20 



APPENDIX 


1 + 


Mg 


2 

m= 


1 


^m,2 

^0,2 



= 0 


(A5) 


If the choice of equations with n =. 1 and n = 2 is made so that S Mg, subtraction 
of equation (A5) from equation (A4) yields 



(A6) 


or 


Mj^-1 

2 (^m+1,1-^0,2 ■ '^m+l,2%l)‘^0!*” = ® 

m=0 


(A7) 


Note that if the mth power on c„ does not exist in the nth equation, then A is zero 
for that particular case. 

Consider equation (A3) first with n = 1 and then with n = 2 so that 


and 



+ 


Mj-1 


I 

m=0 


_jnJ,c “=0 

t a 

^Mi,l 


(A8) 


M, 


Mg-l 


a 


^m,2 

n^O ^Mg,2 


I 


c “= 0 
a 


(A9) 


If Mj = Mg, subtraction of equation (A9) from equation (A8) yields 


Mj-1 

I 

m=0 


^m.l 

^Mj,! 


'^m,2\ m 


= 0 


(AlO) 
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or 


Mj-1 

^ (V,1^M2,2 

m=0 ^ ^ 




A 


m 


a 


= 0 


(All) 


For the case where = Mg, equations (A7) and (All) represent two independent 
equations in which the highest order of c^ is Mj - 1 rather than Mj as it was in 
equation (Al) for n = 1 and n = 2. However, if > Mg, equation (A7) is still 
obtained, but instead of using equation (All), equation (Al) with n = 2 is used unchanged. 


Mg 

I 

m=0 


A „c = 0 
m,2 o? 


(A12) 


For this case where Mj^ > Mg, equations (A7) and (A12) represent two independent equa- 
tions in which the highest order of c^ is Mj - 1. In either case, the coefficients in 
equation (A7) and in equation (All) or (A12) are polynomials which can be redefined for 
convenience so that 


Mj-1 




B ,c"^ = 0 
m,l a 


(A13) 


and 


Mg 

I 

m=0 


B «c ^ = 0 
m,2 a 


(A14) 


The procedure to this point is now repeated by starting with equation (Al) but with 
a different pair, for example, n = 3 and n = 4. If K is odd, one of the n original 
equations indicated by equation (Al) will be used twice so that only one independent equa- 
tion is developed from that pair. A general expression for these reduced equations, 
which includes equations (A13) and (A14) as special cases, can be written as 

M„-l 

I VnS"’=“ (n=l,2,...K) (A15) 

m=0 
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Because this equation is of the same form as equation (Al), it follows that the order of 
these K equations can ultimately be reduced to first order with respect to c^ that is, 

%n^a + Qo,n = ° (n = 1, 2, . . .K) (A16) 

where and Qq ^ are involved polynomial functions of K - 1 mass-fraction 

variables excluding c^. 

The variable c^ can be eliminated from equation (A16) to give the following 
K - 1 equations which do not contain c^ 


%,n^l,n+l 


" ‘^l,n^0,n+l 


= 0 


(n= 1,2,. . .K-1) (A17) 


Equation (A17) can now be written in the form 

I Am,n<=r=“ (n = 1,2,. . .K-1) (A18) 

m=0 


wherein c^ is the mass fraction of one of the remaining K-1 species. 


highest order of c^ for the nth equation, and A^ 




is the 


are polynomial functions of 
K - 2 mass-fraction variables excluding c^, as well as functions of the known quantities 




a, 


i,k’ 




^i,j, 


\,k+l* 


Equation (A18) represents K-1 equations which contain K-1 mass-fraction 
variables. The reduction procedure used to obtain equation (A18) which represents 
K-1 equations and unknowns by starting with equation (Al) which represents K equa- 
tions and unknowns can be used as many times as is necessary to obtain a single poly- 
nomial expression involving only one mass-fraction variable, for example, Cj^, 



(A19) 
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